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ABSTRACT 

We carry out high-resolution adaptive mesh refinement simulations of a cool core cluster, resolving 
the flow from Mpc scales down to pc scales. We do not (yet) include any AGN heating, focusing 
instead on cooling in order to understand how gas gets to the supermassive black hole (SMBH) at 
the center of the cluster. We find that, as the gas cools, the cluster develops a very flat temperature 
profile, undergoing a cooling catastrophe only in the central 10-100 pc of the cluster. Outside of 
this region, the flow is smooth, with no local cooling instabilities, and naturally produces very little 
low-temperature gas (below a few keV), in agreement with observations. The gas cooling in the center 
of the cluster rapidly forms a thin accretion disk. The amount of cold gas produced at the very 
center grows rapidly until a reasonable estimate of the resulting AGN heating rate (assuming even 
a moderate accretion efficiency) would overwhelm cooling. We argue that this naturally produces a 
thermostat which links the cooling of gas out to 100 kpc with the cold gas accretion in the central 100 
pc, potentially closing the loop between cooling and heating. Isotropic heat conduction does not affect 
the result significantly, but we show that including the potential well of the brightest cluster galaxy 
is necessary to obtain the correct result. Also, we found that the outcome is sensitive to resolution, 
requiring very high mass resolution to correctly reproduce the small transition radius. 



1. INTRODUCTION 

The majority of baryonic matter in galaxy clusters re- 
sides in the for m of virialized h ot gas that emits in the 
X-ray band (e.g. Lin et al.|2003 ). Observations of the X- 
ray surface brightness, along with imaging spectroscopy, 
allow the measurement of the density and temperature 
of the intracluster medium (ICM). Many, if not most, 
relaxed clusters possess a cool core (CC), where the tem- 
perature drops by a factor of 2 to 3 compared to the 
cluster outskirts and the radiativ e cooling time is much 
shorter than the Hubb le time (e.g. Cavagnolo et al.|2008 
Sanderson et al. 2009 ) . In a steady state, a "cooling flow 



of 100s M f7)/yr would de velop in a typical CC cluster (see 
review by Fabian 1994). However, XMM grating spec- 
tra have revealed a dearth of gas at temperatures lower 
than about one-third of the cluster viria l temperature, 
despite the short cooling time of the gas ( Tamura et al. 
2001[|Peterson et al.|2003"l|Sanders et al.|2008[). In addi- 
tion, the observed star-formation rate (e.g. Cardiel et al. 
1998 1 is typically a factor of ~ 10 — 100 times lower than 
the mass deposition rate predicted by classic cooling flow 
models. This discrepancy implies that there is some pro- 
cess that heats the gas, preventing it from cooling. Many 
mechanisms have been proposed, includi ng thermal con- 
duction (e.g. Zakamska & Narayan 2003 ) , active galactic 
nuclei (AGN) feedback (e.g. Bohringcr et al. 1993; Mc- 
Namara et al .||2005| |BrighentT&: Mathew s 2006), MTTLT 
d cosmic rays combined with c onduction 

20081), and tur- 



wave-mediate* 
(e-g 



Loewenstein et al. 1991 Guo & Oh 



bulence co mbined with conduction (e.g 



ennis & Chan- 



dran|2005 |. 

AGIN feedback is considered to be the most plausible 
heating sour ce for a number of reasons (McNamara & 
Nulsen|2007). First, there is an observed correlation be 



tween the presence of cool cores and signs of AGN ac- 
tivity (e.g. Dunn fc Fabian|2 006). In addition, in nearby 
CC clusters, X-ray cavities or bubbles are observed in the 



core (e.g. Fabian et al. 2000) and are thought to be in- 
flated by the interactions between powerful jets and the 
surrounding gas. Finally, the estimated energy injected 
into the surrounding gas by the bub bles can balance cool - 
ing in most of these systems 



(e.g. 



Birzan et al. 



2004) 
ler the 



However, how the gas feeds the AGN, i.e. whet 
gas accretion is dominated by the "h ot" mode described 
in the Bond i-Hoyle accretion model (Bondi 1952 Allen 



et al. 2006) or by the "cool" mode assuming that the 
hot gas first fragments within the sonic radius and falls 
onto the black hole as cold clouds (e.g. |Balbus fc Soker 
1989), and how AGN deposits the energy to the ICM 
(whether kinetic or t hermal feedback dom inates) are still 



open questions - see Croton et al. ( |2006 ) and references 



therein. To address these issues, we need to first under 
stand exactly how gas cools and eventually accretes onto 
the black hole. 

Theoretical work has been carried out using linear per- 
turb at ion_thej3ry_toJri^^ 

ity (e.g. |Binney et al.|2009||Balbus fc Soker|1989[ ), deriv- 
ing the gl obal dynamics of the coolin g flow in a steady 
state (e.g. Quataert fc Narayan | 200 |) , or simulating the 
global cooling how numerically (e.g.|Cattaneo fc Teyssier 



2007| |Gaspari et al.||2010] |Guo fc Mathews||2010[ ). How- 



ever, due to its complex and non-linear nature, an accu- 
rate picture of the cooling flow cannot be derived from 
analytical models, and the previous simulations usually 
suffer from a low resolution with finest grid sizes typically 
around a kpc. 

In this paper, we carry out 3D simulations using the 
adaptive-mesh hydro code Enzo to investigate the onset 
of the cooling flow. Our primary objective is to explore 
the cooling catastrophe in detail and examine how gas 
cools and flows inwards in the cluster center. We do 
not attempt to model the accretion on to the SMBH it- 
self and the resulting feedback, and so our model is not 
a complete description of thermally-balanced feedback; 
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however, by examining the detailed evolution of a clus- 
ter as it undergoes a cooling catastrophe, we can get a 
better idea of the evolution and production of cold gas 
in the cluster core. Key issues that this study aims to 
address include: (i) confirmation of linear calculations 
that local cooling instabilities do not grow significantly 
before a cooling catastrophe occurs, (ii) the location and 
amount of cold gas produced by the global thermal in- 
stability, (iii) the rate of gas accretion on to a central 
SMBH, (iv) observational signatures during the cooling 
catastrophe (in particular, the lack of cool gas observed 
in X-rays), (v) the impact of other processes (such as 
isotropic thermal conduction and Type la SN heating 
rates) on the cooling catastrophe. 

This paper is organized as follows, in Section [51 we 
describe our code, initial conditions and (non-standard) 
refinement method; in Section [31 we present the main 
results from our standard simulation; in Section [4j we 
discuss the related physical and numerical issues, and 
compare our results with observations and previous work. 
We conclude in Section [5) 

2. METHODOLOGY 

The simulations are performed using the Enzo Adap- 
tive Mesh Refinement (AMR) code, a parallel, Eulcrian 
hydrodynamics scheme. An important advantage of us- 
ing an AMR code is that it achieves high levels of resolu- 
tion by only refining areas that need it and thus reduces 
computational time. The code initially places a uniform 
root or "parent" grid consisting of relatively large cells 
over the whole simulation box. If further resolution is 
required, a finer "child" grid is placed inside the parent 
cell and the properties of each of its grid cells are then 
computed. The child grid then becomes a parent grid at 
the next step and this process can be repeated until the 
desired level of resolution is reached. A more detailed 
descrip tion of Enzo can be fou nd in |Bryan fc Norman| 
(1997) and |Q'Shea et al.| (2004), and references therein. 

The hydrodynamics algorithm we choose for most of 
our runs is a three-dimensional version of the Zeus as- 



trophysical code developed by Stone & Norman ( 1992 1 
in Cartesian coordinate system. It is a simple, fast, and 
robust algorithm that allows large problems to be run 
at high resolution. For comparison, we also perform 
test runs using the piecewi se parabolic method (PPM) 
( [Woodward fe Colella||l984] ) . 

In the following, w e de scribe the initial setup of the 
sim ulatio ns in Section |2.1| our refinement strategy in Sec- 
tion |2.2[ and the physics included in the simulations in 



Section 1273 

2.1. Initial Conditions 

Our simulations start with an isolated idealized galaxy 
cluster roughly in hydrostatic equilibrium. We adopt the 
Perseus Cluster, a massive relaxed CC cluster that has 
been well observed, as our model template. Initially 
the gas is taken to be spherically symmetric and the 
azimuthally averaged temperature profile is taken from 
analytic fits to the ob servations of the Perseus Cluster 
( |Churazov eraI1[2004| : 



T : 



1 + (r/71) 3 
2.3 + (r/71) : 



keV 



(1) 




10" 



_t i i i ■ i 



0.01 



0.10 



1.00 10.00 100.00 1000.0C 

r (kpc) 



Fig. 1. — The gravitational acceleration for our idealized cluster 
(solid line) along with the contribution from different components. 
All contributions except for the gas are static and do not change 
during the simulation. 



where r is the radius to the center of the cluster in kpc. 
This temperature profile is flat for r > 0.5 Mpc, while 
Perseus, like ma ny clusters, has a decli ning profile at 



large radius (e.g. Vikhlinin et al. ||2006[ ). Since we are 
interested in the evolution in the center (r < 0.2 Mpc), 
we do not expect this difference to be signi ficant. The 
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The power-l aw index of the last term is slightly steepened 
compared to Mathews et al. ( 2006 ) so that at small radii 
(up to a few hundred kpc) the gas density profile does not 
significantly deviate from that obtained from the X-ray 
observation, but at large radii, where the fit is poorly 
constrained, the density drops as r~ 2 ' 2 , which is more 
consistent with cosmological simulati ons and fits better 



with the NF W dark matter profile (Navarro, Frenk & 
White 19961. We compute the initial pressure from the 



density and temperature assuming the ideal gas law with 
7-5/3. 

The gravitational potential is the combination of the 
self-gravity of the ICM (which does not dominate but is 
included for completeness), and three static components: 
the dark matter halo, the stellar mass of the brightest 
cluster galaxy (BCG) and the super-massive black hole 
(SMBH) in its center. The dark matter halo follows a 
NFW distribution that matches the observed gas density 
and temperature of Perseus: 



NFW 



(r) = 



Po 



(f-)(l + FT) 5 



(3) 
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where pg FW = 7.5 x 10 14 M /Mpc 3 , and r s = 0.494 
Mpc is the scale radius. Since the gas density and tem- 
perature are only observed up to a few hundred kpc and 
are uncertain at large radii, the NFW parameters can 
vary depending on the outer radius we choose when fit- 
ting the model. We experimented with slightly different 
NFW parameters and the results were not affected. 
The stellar mass profile of the BCG is: 



G 



p 0.5975 



3.206 x 10- 



„1.849 



1.861 x 10- 



in_cgs units with s = 0.9 and r in kpc (Mathews et al. 
20061. The SMBH is treated as a point mass at the 



-l/s 



Ml 



very center of the cl uster with Msmbh = 3.4 x 10 M Q 
(|\\ ilmaii <n aTffioi} . 

The self-gravity of the ICM is computed at each step, 
but it does not have much impact in our major runs be- 
cause the total mass of the cool gas in the cluster center is 
small compared to the BCG and SMBH. The mass of the 
accumulated gas in the center of the core only becomes 
significant at late time in our low resolution sim ulati ons 
where the cooling flow is stronger (see Section 4.5 for 
discussion) . The gravitational acceleration from different 
components is plotted in Figure [l] As can be seen from 
that figure, the radius of influence of the SMBH (where 
its acceleration dominates) occurs at approximately 50 
pc. 

In order to give the gas some initial angular momen- 
tum, we set the cluster initially rotating slowly around 
the z axis in our simulation. The angular momentum of 
a galaxy cluster can be characterized by its spin pa ram- 



Peebles|(fl97l1) and 



eter A, which was first introduced by 
later interpreted also as A = oj/lu sup (Gottlober &; Yepes 
2007 1 , where ui sup is the angular velocity that would pro 
vide rotational support to the system. As a simplifica- 
tion, we use the same rotational velocity v at all radii and 
to v/r s . Thus the spin parameter A can be expressed 

as A = v/v c , where v c = \J G ^ Is with M s being the to- 
tal mass inside the scale radius r s . Typical values of A 



for a cluster are around 0.05 (e.g. Bullock et al. 2001 ). 
For our model cluster, we take a constant rotational ve- 
locity of v = 50 km/s, which is consistent with cluster 
simulations. 

We also give each cell an initial random velocity that 
could potentially seed small-scale instabilities. The ran- 
dom velocity obeys a Gaussian distribution with a stan- 
dard deviation of 200 km/s for each of the three compo- 
nents, a typical value for the turbulent m otion of gas in 



galaxy clusters ( Ruszkowski fc OhpOlO ) 

To set up the initial cluster configuration, we employ 
an iterative technique: starting with the root grid, we 
use the equations above to set the density, temperature 
and velocities. We then apply the refinement criteria 
discussed in the next section to obtain additional levels 
of refinement, applying the initial conditions to each level 
and reapplying the initial conditions. The result is a grid 
hierarchy which is self-consistcntly refined, and contains 
high-resolution initial data. 

2.2. Refinement Strategy 

We refine a cell whenever any of the following three 
criteria are met: 



(i) The gas mass criterion - a cell is refined if the gas 
mass in any cell exceeds 0.2 times the gas mass in one 
cell of the root grid, or more precisely, when: 



m ce ll > 0.2 p mC an(^ ) 3 X 2 

JVroot 



al 



(5) 



where p mea n is the mean density of the universe at the 
initial time, L = 16 Mpc is the co-moving box size, I is 
the level of refinement, and N root is the number of cells 
on the root grid in each dimension. We use N IO ot = 256 
for our standard simulation. This defines a critical gas 
mass that is always refined, which we then modify (with 
the factor 2 al ) depending on the level of refinement. If 
a = 0, the refinement is Lagrangian and the limit for 
TO C eii does not change with refinement level in this case, 
a given mass clump will be resolved by the same number 
of cells as it collapses to smaller sizes. A negative value 
of a makes the refinement "super-Lagrangian" and the 
maximum value of m ce u decreases with I. The decrease is 
faster as a becomes more negative. In our simulation, we 
find that choosing a more negative a can better res olve 
the early evolution of the cooling flow (see Section 4.5 
for discussion). We use a = —1.2 for our standard run, 
which reduces the maximum mass of a cell by a factor of 
3.8 x 10 -6 going from I = to I = 15. In the standard 
run, the level of refinement in the initial conditions is 
/ = 8, which gives a maximum cell mass (i.e. roughly the 
mass resolution in the initial conditions) of about 8 x 10 5 

Mq. 

(ii) The cooling criterion - we also apply refinement 
whenever the ratio of cooling time to sound-crossing time 
over the cell becomes too small (i.e. t coo \/t cross (Ax) < 
(3). The isobaric cooling time for an optically-thin plasma 
is given by 



t, 



cool 



|nfcfcT 
n 2 A(T) ' 



(6) 



where kb is the Boltzmann's constant, n — p/ pmn is the 
particle number den sity and A(T) is the cooling function 
( Schure et al.| [2009). The sound crossing time over a cell 
is icrossi^^) = Ax/c s , where Ax is the cell size, and 
c s = y/jP/p is the sound speed. We use this criterion 
because when the cooling time is shorter than the sound- 
crossing time, the gas drops out of pressure equilibrium, 
and the cooling becomes isochoric. We use the limit (3 — 
6, which is a somewhat arbitrarily chosen value larger 
than 1 so that we can fully resolve this transition. We 
have experimented with higher value of (3 = 12 and (3 — 
20 for this ratio and the results do not change. 

(iii) The Jeans length criterion - a cell is refined when- 
ever its size is larger than 1/4 of the Jeans length, to 
ensure that we properly reso lve any gravitational insta- 



b- 
ilities (|Truelove et al.|r997|) 



In our simulations, the baryon mass criterion is impor- 
tant at early stages, well before the cooling catastrophe 
occurs, while the cooling criterion becomes important 
later as the gas density grows in the center, and the 
radiative cooling becomes more significant. The Jeans 
criterion is never important for our standard runs. 

2.3. The Simulations 

All simulations in this paper include radiative cooling. 
The radiative cooling function we use is adapted from 
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Fig. 2. — The evolution of gas density, temperature and pressure. In the top row, these are: (left) spherically averaged gas density 
profiles, and (right) the time evolution of the central gas density. On the bottom row, we plot (left) the mass weighted temperature profile, 
and (right) the gas pressure profile. The black dashed lines correspond to the initial conditions. Each solid line is from one output, with 
different colors corresponding to the different times on the right panel. 



the SPEX package (Schure et al. 2009), assuming equi- 
librium cooling. Non-equilibrium effects will play a role 
in cooler gas, but as this only occurs in the very center 
of our simulations, we do not expect this to play an im- 
portant role. We a dopt a metallicity of o ne-half solar for 
our model cluster ( Schmidt et al. |2002 1 . The cooling is 
truncated at temperatures lower than 10 4 K , because we 
are not interested in following the evolution of gas below 
this temperature. 

Our standard run, which will be the primary simula- 
tion analyzed in this paper, has N mot = 256, with a 
co-moving box size of L = 16 Mpc, and a maximum 
refinement level of Z max = 15. The minimum phys- 
ical size of a grid cell in this simulation is given by 
Acc m i„ = L/(A r root 2' max ), which is about 2.0 pc. 



section 
We a 



To test the dependence of the results on the resolution, 
we also performed simulations with N root = 64 , 128 and 
256. We also carry out simulations with the same N root 
but diff eren t values of a. These runs will be discussed in 
4.5| See Table [T] for the parameters of all runs, 
so performed several runs with isotropic thermal 
conduction suppressed from the standard Spitzer value 
by different factors of / con d in order to test its influence 
on the evolution of the cooling catastrophe (see Section 
4.3 for details). Two different suppression values were 
tested (/cond = 0.1 and / con d = 0.3). We did not in- 
clude star formation or any feedback mechanism in our 
simulations. 

In order to facilitate analysis of the simulation, Enzo 
generates an output when a new refinement level is 
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TABLE 1 
Simulation Parameters 



Simulations N r , 



Figures 



1 

2 

3 (standard) 

4 

5 

6 

7 



64 

128 

256 

128 

128 

128 

128 




reached for the first time. This produces outputs 
throughout the collapsing phase up until the point when 
the maximum refinement level l max is achieved. To bet- 
ter resolve the evolution after this, starting from a few 
Myr before l m ax is reached, the simulation writes one 
output every 1/3 Myr and continues to do so afterwards. 
We run the simulation until the estimated energy from 
AGN feedback is more than strong enough to offset cool- 
ing, which occurs approximately 10-20 Myr after reach- 
ing ^ max - 

For completeness, we note that these simulations were 
run using comoving coordinates, but over a sufficiently 
brief span (Az < 0.2) that expansion did not play a role. 
We adopted cosmological parameters corresponding to 
Ho = 50 km s _1 Mpc -1 , fl\ = and O m = 1, but note 
that the cosmological parameters are irrelevant to the 
physics inside the self-gravitating cluster and thus have 
no significant influence on its evolution or on the initial 
condition (they just serve to set the internal code units). 
The simulation starts at z — 1, which is chosen to give 
enough time for the cluster to evolve. 

3. RESULTS 

In this section we present and discuss the results from 



our standard simulation. In Section 3.1 
time evolution of the cooling gas. In Section |3.2 



we describe the 
we 



describe when and how the cooling catastrophe happens. 

3.1. Cluster Evolution 

We show in Figure [2] the evolution of the gas density of 
our standard simulation with N root = 256. Throughout 
the cluster core (r < 30 kpc) the density slowly increases 
and the temperature decreases over the first few hun- 
dred million years. The density in the center (r < 100 
pc) grows more compared to the outer part of the core 
region r ~ 0.3 — 10 kpc, where the temperature shows a 
plateau. After about 300 Myr, which corresponds to the 
cooling time of the gas in the initial conditions, the evo- 
lution rapidly accelerates, and the density increases by 
several orders of magnitude within a million years. No- 
tice that the outputs are not evenly spaced in time (out- 
put times are marked and color-coded in the upper-right 
panel of Figure 2]) This marks the onset of the cooling 
catastrophe, andoccurs in the central region of the clus- 
ter. We refer to the outer boundary of the region where 
this catastrophe happens as the "transition radius" . This 
is different from the traditional cooling radius which is 
defined as the radius at which t coo i is equal to the Hubble 
time (and which occurs near 100 kpc). The huge jump 
of the central density from ~ 10~ 23 g cm -3 to ~ 10 -21 g 
cm -3 happens within only a short period of time, about 
3 Myr, which also corresponds well with the drastic de- 
crease in the central temperature seen in the lower-left 
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Fig. 3. — Radial velocity of the gas (spherically averaged), color- 
coded for the same times as in Figure [ 



panel of Figure [2] After the cooling catastrophe, the 
transition radius steadily grows in time, although points 
exterior to the transition radius evolve only very slowly. 

The pressure profile (Figure [2j lower right panel) shows 
relatively little variation over time, at least in the outer 
parts. As more gas flows to the cluster center, the pres- 
sure slowly builds up, although the slope in the interme- 
diate region (between 100 pc and 10 kpc) is quite shallow. 
As the cooling catastrophe occurs, the pressure in the 
very center rises quickly. A small region with a slightly 
positive pressure gradient forms later right at the inner 
edge of the transition radius (at r ~ 50 pc), creating a 
slight pressure "hole" (or inversion). This pressure hole 
does act to drive some gas inflow, however, the pressure 
profile is remarkably constant at the transition radius, 
while the density and temperature change by orders of 
magnitude. Both the size and the depth of the pressure 
hole are sensitive to the resolution of the simulation, as 



we will discuss in Section 4.5 becoming less prominent as 
the resolution increases. After the cooling catastrophe, 
the pressure profile is nearly time-independent. 

It is clear from Figure [2] that the cluster can be di- 
vided into three regimes in radius. In the outer region 
(r > 20 kpc), where the NFW dark matter dominates 
the gravitational potential, the gas properties stay rel- 
atively constant during the simulation. Inside the core, 
but outside the transition radius (0.1 kpc < r < 20 kpc), 
where the gravitational potential is dominated by the 
BCG, the density and pressure increase slowly while the 
temperature profile is nearly flat. Finally, in the very 
center of the cluster (r < 0.1 kpc), within the radius of 
influence of the SMBH where the gravity is dominated 
by the SMBH (~ 50 pc, see Figure RU, the cooling catas- 
trophe first occurs at small radius, then the transition 
radius moves slowly outwards with time, from ~ 10 pc 
when the collapse first happens to about ~ 100 pc, about 
10 Myr later. 

Figure [3] shows the evolution of the radial inflow ve- 
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Fig. 4. — The gas surface density (upper left), X-ray weighted temperature (upper right), pressure (lower left) and the X-ray luminosity 
(0.5 — 2 keV) (lower right) in the cluster core projected along the direction of the initial angular momentum z-axis at late stage of the 
simulation, about 296 Myr after the initial time. The size of the region is roughly (16.6 kpc) 3 . 
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locity of the gas. Initially it is zero, and in the ab- 
sence of cooling stays that way except in the outer re- 
gion (r > 100 kpc), where the cluster is not in perfect 
hydrostatic equilibrium^] As radiative cooling proceeds, 
a slow but steady inflow develops. This is slow (~ 20 
km/s) and nearly constant over a large range in radius, 
from 20 kpc to a few hundred pc. Notice that this is 
not a steady state cooling flow because a constant veloc- 
ity implies an accretion rate that rises with radius - we 
discuss this in more detail below. As the cooling catas- 
trophe takes hold (recall that the outputs are not evenly 
spaced in time), the inflow velocity in the central region 
rises up to a few hundred km/s inside 100 pc, and then 
drops after an accretion disk forms. 

So far, we have focused on one-dimensional analysis 
of the cluster evolution. In fact, the evolution is remark- 
ably symmetrical, despite the fact that we seed the initial 
conditions with substantial small-scale velocity motions 
which quickly generate density and temperature pertur- 
bations. We show projections of the density, (X-ray 
weighted) temperature, pressure and X-ray luminosity 
of the central 16.6 kpc of the cluster in Figure [4] These 
are shown 296 Myr after the initial time, shortly (~ 16 
Myr) after the cooling catastrophe occurred. Clearly the 
initial seeds have damped and the profiles are remarkably 
uniform. A close examination of the temperature image 
reveals faint spiral structures due to the slight rotation 
imprinted in the initial conditions. This demonstrates 
that the flow does not exhibit local thermal instabilities, 
as predicted by linear perturb ation theory (Malagoli et 
1987] |Balbus & Soker||1989[ ) 



10.0 



al 



he primary instability (or catastrophe) which grows 
is a global one, as discussed above. This can be seen 
as the increase in density in the very center of the pro- 
jected image (and the decrease in the temperature). The 
pressure map is quite smooth. The X-ray map in Fig- 
ure [4] shows a small rise in the central emissivity in this 
region, but it can be appreciated that this would be dif- 
ficult to observe. We return to t he o bservability of this 
lower temperature gas in Section [4. 2| 

We stop the simulation about 20 Myr after the cool- 
ing catastrophe first occurs because the estimated energy 
output from AGN feedback rapidly exceeds the energy 
loss through radiation. We estimate that the feedback 
energy from the AGN as 



Eagn — eMiQQC 



(7) 



where M400 is the mass accretion rate measured at r = 
400 pc, which is just outside the region dominated by 
the potential well of the SMBH, and e is an efficiency 
parameter discussed in more detail below. We stop the 
simulation shortly after this feedback rate exceeds ixsoj 
which is the total calculated X-ray luminosity from the 
cooling gas inside r < 50 kpc. 

The parameter e is the efficiency relating this mass ac- 
cretion rate to the feedback energy that heats the ICM, 
and is really the product of three somewhat poorly con- 
strained efficiencies. The first is the fraction of M400 
which actually accretes on to the black hole. Because we 

1 This occurs because we use observationally-motivated density 
and temperature pro files which are not in perfect hydrostatic equi- 
librium. See Section 14.61 for a discussion. 
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Fig. 5. — The estimated ICM heating rate from the AGN (com- 
puted with Equation[7] shown as the solid line), compared with the 
total X-ray luminosity in the central 50 kpc region (dashed line). 



measure this value quite close to the SMBH-dominated 
region (at 400 pc), we argue that a significant fraction 
(at least 10%) of this mass will accrete on to the black- 
hole. The second unknown efficiency is that associated 
with the conversion of accreted mass into energy, which 
we take to be approximately 10%. Finally, because the 
accretion rate is likely to be sub-Eddington, this feed- 
back will be so-called 'radio-mode', and we assume that 
a large fraction of the feedback energy ends up heating 
the ICM. For the product of these three efficiencies, we 
take a conservative value of e = 0.01, although we recog- 
nize that the final value is uncertain. 

As shown in Figure [5j as the gas inflow increases, 
Eagn increases and exceeds Lxso after the cooling catas- 
trophe. At the last step, Eagn ~ 2.6 x 10 44 ergs/s and 
Lx50 ~ 1.5 x 10 44 ergs/s, which has increased only about 
50% compared to the initial L X 5o- We stop our simula- 
tions at that point because AGN feedback, which is not 
included, would be important and we only want to focus 
on the onset of the cooling flow. The total X-ray lumi- 
nosity of the cluster is 6.8 x 10 44 ergs/s at the end of our 
simulation. 

3.2. The Global Cooling Catastrophe 

To better understand the physics behind the cluster 
evolution, we examine several relevant timescales in Fig- 
ure [6j The dynamical time of a gaseous system can be 
defined as: 

2 ^/GM dyn (r) 

where Md yn (r) is the enclosed mass of the system. 
The system undergoes gravitational collapse when £dyn 
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is shorter than the sound crossing time 

v 

£sound(V) = — , (9) 

where r is the radius of the region. 

Figure [6] shows the cooling, dynamical and sound cross- 
ing time scales before (left panel) and after (right panel) 
the cooling catastrophe happens. Initially the cluster is 
stable with t coo \ > t^ yn > i SO und- As the gas condenses 
and cools nearly isobarically in the center, the cooling 
time t coo i decreases while i SOU nd slowly increases (both 
at fixed radius). When i coo i drops below td yn , which first 
occurs around r < 50 pc, the gas undergoes a full-fledged 
cooling catastrophe. At the same time, i S ound exceeds 
£dym making the core region gravitationally unstable and 
causing the gas to collapse to the cluster center under 
gravity. Meanwhile, the gas can no longer maintain pres- 
sure equilibrium since t coo \ < t SO und- Also shown in Fig- 
ure [6] are the thermal conduction and the Supernovae 
la (SNI a) h eating timescales which will be discussed in 
Section |4~3| and Section |4~4{ respectively. 

This analysis leaves open two questions. First, why 
does the gas temperature remain nearly constant in the 
intermediate region, between the transition radius (at a 
hundred pc) and the start of the cool core (at r ~ 20 
kpc)? Second, what supports the gas inside the transi- 
tion radius? 

The first question arises because of the presence of the 
large temperature plateau noted earlier, which persists 
despite the fact that we evolve the system for many cool- 
ing times at these radii. Clearly, the cooling must be bal- 
anced by some form of heating, and since the gas motions 
are subsonic, the only real candidate is compression heat- 
ing. We estimate the compressional heating timescale as: 

icomprcssM = ^ _ _ " ) (10) 

where 7 = 5/3 and V • v = i d ^ r g ^ ■ For simplicity, we 
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Fig. 7. — The ratio of the (radial) compressional heating 
timescale over cooling time; see text for definitions. 

only include the radial term, focusing on the role that 
the cooling flow plays in establishing the temperature 
plateau (we note that non-radial terms could decrease the 
heating time in the center, although manual inspection 
of the velocity indicates that the effect is minor). We 
plot in Figure [7] the ratio of tcompress over t coo \ . At early 
times, over a large range of radius, t CO mpress ~ icool- The 
inflow velocities which drive this heating are slight, a few 
tens of km/s over much of this range, as seen in Figure [31 
rising at small radius as the global cooling catastrophe 
sets in. 
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Fig. 8. — The circular velocity of the gas (dashed lines) compared 
to the Keplerian velocity (solid line) . The gas becomes rotationally 
supported in the very center at late times. 



This brings us to the second question, the evolution 
of the gas within the transition radius. As the cluster 
evolves, the cooling rate becomes higher in the center, 
requiring a larger inflow velocity to provide sufficient 
compressional heating. When this required velocity ex- 
ceeds the sound speed, or if the gas becomes rotationally 
supported, compressional heating can on longer balance 
cooling and i C ompress becomes larger than i coo i inside the 
transition radius. If the inflow velocity grew to the sound 
speed, as would occur for a purely radial evolution, the 
transition radius would be identified as a sonic point, and 
inside the gas would freely fall toward the black hole. 
We found this to occur in our lower resolution simula- 
tions (discussed in more detail below); however, in our 
best resolved, standard simulation, we find that the cold 
gas inside the transition radius forms a rotationally sup- 
ported disk. This is shown in Figure [8j which shows an 
estimate of the rotational velocity of the gas (computed 
by dividing the magnitude of the total specific angular 
momentum of the gas in a shell by the radius of that 
shell), compared to the Keplerian velocity. At late times, 
inside the transition radius, the gas becomes rotationally 
supported. 

In Figure|9j we show the projected density and density- 
weighted temperature in the central 330 pc for a slice of 
gas with a z-thickness of about 16.6 pc. This z-projection 
clearly shows the disk (x- and y-projections - not shown 
here - clearly demonstrate that this is a thin disk), which 
has a radius of about 50 pc. In fact, in this particu- 
lar run, we find an inner disk and an outer polar ring, 
which shows some sign of forming denser fragments; typ- 
ical densities in the disk are of order 10 3 cm -3 , but note 
that the disk is not well-resolved and so the detailed disk 
structure shown here should be treated with caution. 



4. DISCUSSION 

In this paper, we have presented results from the high- 
est resolution simulation of the onset of cooling in a cool- 
core cluster. We have demonstrated that the flow is re- 
markably uniform, with thermal instabilities not growing 
outside of the central few hundred pc, where the tem- 
perature drops rapidly at a point we have termed the 
transition radius. In addition, we have shown that the 
flow natufrally generates a nearly constant temperature 
state outside of this transition radius. Inside, a rotation- 
ally supported accretion disk forms around the central 
SMBH. This time-dependent flow is not in steady state 
and is not, without additional heating, a solution to the 
cool core problem. Nevertheless, we have made substan- 
tial progress in delineating exactly when and where cold, 
dense gas first condenses out of the flow. 

However, there are a number of unanswered questions, 
including a bet ter understanding of why this solution 
occurs (section 



4.1), a detailed examination of the ob- 



servatio nal p redictions of the the final simulation state 
(section |4.2| , a first attemp t to examine the impact of 
thermal conduction (section 4.3) a nd T ype la SN h eat 



4.5 



to 



ing from stars in the BCG (section 4.4). In section 
we show that high numerical resolution is required 
obtain these results, and with lower resolution, the tran- 
sition radius first forms at much larger radius. Finally, 
we argue that the results are robust to changes in the 
initial conditions (section 4.6), a nd t hen compare these 
results to previous work (section 4.7). 



4.1. Cluster Profile 

To better understand what determines the structure of 
the gas and why there are three regimes seen in the gas 
density, temperature and pressure profiles in Figure [2j 
we carry out an approximate analytic analysis assuming 
hydrostatic equilibrium, which is valid before the cooling 
catastrophe happens (or at radii larger than the transi- 
tion radius) when the inflow velocity v r -C c s . 

In hydrostatic equilibrium, the gas pressure P(r) bal- 
ances the gravitational acceleration g(r) 



dP(r) 
dr 



-p{r)g{r) 



(11) 



where P = pkbT / pmn and p(r) is the gas density. This 
by itself is not sufficient to determine the density and 
temperature profiles uniquely; however, if we also use 
the fact, as suggested by Figure [7J that compressional 
heating balances cooling, we can make more progress. 
We write this as, 



pgv r 



2 A{T) 



(12) 



where A(T) is the cooling rate and v r is the radial veloc- 
ity. We have assumed the work done on a fluid element 
as it flows inward is balanced entirely by radiative cool- 
ing. Simplifying further, if we assume that A(T) » Ao 
is independent of temperature, as is appropriate for gas 
around a few keV, in the cool region (see Figure [2]), and 
that the inflow velocity is nearly constant (see Figure [3]), 
then we can obtain an expression for the gas density: 
n = pmiigv/ Ao, and so p oc g. 

We expect this result to hold primarily in the inter- 
mediate region of the flow, from the transition radius at 
about 100 pc to about 20 kpc. This region is dominated 
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Fig. 9. — The projected density (left) and density-weighted temperature (right) of a slice of gas in the very center of the cluster at the 
same time as in Figure [4] The projection is along z-axis and the size of the region is about 330 pc x330 pc X16.6 pc. 



by the potential of the BCG, which is reasonably well 
fit over this range with a power law g(r) ~ r o"°_ 
predicted by the above argument, this is an excellent ap- 
proximation to the density profile seen in Figure [2j If we 
combine this result with the equation of hydrostatic equi- 
librium, we predict that the temperature profile should 



be T(r 



which is also in reasonable agreement 



with the nearly flat temperature profile we find in Fig- 
ure [2] over this radial range. 

To see how changing the form of the gravitational po- 
tential in the cluster core can affect the structure of the 
cooling gas, we performed a test simulation with only 
the NFW dark matter but without the BCG and the 
SMBH and compare it with the standard run with the 
same number of cells on the root grid (N root — 128). The 
cooling catastrophe starts inside r ~ 1 kpc, much larger 
than the standard run. This is because the tempera- 
ture decreases as T ~ r towards the center throughout 
the core under the NFW gravitational potential and no 
temperature plateau forms. 

4.2. Observational Comparison 

One problem with the classic cooling flow model is 
that the expected amount of cool gas (i.e. gas with 
T < T vir /3, or T < 2 — 3 keV for massive clusters) is 
much higher than that observed in X-ray observations of 
cool core clusters. In our simulations, due to the small 
size of the region where the cooling catastrophe actually 
happens, the amount of cool gas is small and the X-ray 
luminosity does not change much in the core (see Fig- 
ure pi , which would explain the lack of cool gas observed 
in trie X-ray. This can also be seen from the simulated 
X-ray luminosity (Figure H). Note that feedback is still 
required since this will notTast without energy input: the 
amount of cool gas will grow with time and approaches 



the classic cooling flow p redi ction if we let the simulation 
run further (see Section 4.7). 

Since the observed temperature profile of galaxy clus- 
ters is usually obtained by fitting the X-ray spectra as- 
suming a single temperature for the gas, it is hard for 
us to directly compare our temperature profile to those 
observed (since along a given line of sight, t here exists 
gas with a range of temperatures). However, Fabian et 



al. (2006) fitted the X-ray spectra of a long Chandra 



observation of Perseus with a multi-temperature model, 
allowing them to compute the amount of mass in each 
temperature range. We compare the mass distri but ion 
from our simulation with their results in Figure |10[ as 
well as with the prediction from the classic cooling now 
model. We find that the steep drop-off in mass between 
1 — 2 keV in our simulation is in good agreement with the 
observations and is much lower than the classic cooling 
flow prediction. 

Note that the recovery at ~ 0.5 keV in the observed 
mass distribution is due to the observed filamen ts (this 
also probabl y contributes to the 1 keV bin, see |Fabian| 



et al.| (20061), which may be caused by the interaction 
between the cooling gas and AGN feedback that is not 
included in our simulations. The heating from AGN can 
potentially result in an increased amount of gas at ~ 0.5 
keV, and we will study the impact of AGN heating in 
future work. Additionally, there is an excess at around 3 
keV in our simulations and a lack of gas at around 5 keV 
compared to observations. This is likely due to the initial 
condition we choose in our simulations: we start from 
a cool core configuration that matches Perseus' current 
configuration. We suspect that if we started from slightly 
different conditions - say an initial temperature plateau 
that was marginally higher, the fit could be improved. In 
any case, we are not trying to fit Perseus in detail, but 
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Fig. 10. — Distribution of mass at fixed temperatures bins within 
r < 32 kpc. The black solid line shows the result from our simula- 
tion (296 Myr after the start, 16 Myr after the cooling catastrophe); 
the blue dots are the observed distribution within the innermost 1.5 
arcmin radius (~ 32 kpc) and the red dash ed line is the the cla ssic 
cooling flow model prediction, taken from|Fabian et al.|(2006). 
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Fig. 11. — Comparison of the temperature profiles of runs with 
different suppression factors / con d of the heat conduction coeffi- 
cient. Each profile is shown about 2 Myr after the cooling catastro- 
phe first occurs (although the cooling catastrophe occurs at some- 
what different times in each simulation). The black dash-dotted 
line is the initial temperature profile (taken from observations of 
Perseus and extrapolated to the center). 



simply to make the point that the simulation naturally 
produces a dearth of low-temperature (T < 2 keV) gas. 

4.3. Thermal Conduction 

Thermal conduction can potentially delay or even sup- 
press cooling instabilities. However, the existence of 
a magnetic field can decrease the efficiency of thermal 
condu ction below the classical Spitzer value ks ( Spitzer 
1956 1. This suppression can be parameterized by a factor 
of /cond ranging from 0.1 ~ 1 dep ending on the orienta- 
tion a nd topology of the field (e.g. Zakamska & Narayan 
2003 1 . The characteristic timescale associated with ther- 
mal conduction can be estimated as 



tr 



i(r) 



7 



Pr 2 



1 - 1 /cond^sT 



(13) 



7/2 



where ks = 1-84 x 10~ 5 T 5 / 2 /lnA ergs s _1 cm _1 K" 
with the usual Coulomb logarithm In A w 37. 

Thermal conduction is more efficient along the mag- 
netic fields than perpendicular to them. In the center 
of cool core galaxy clusters, where the temperature is in- 
creasing with radius, gas with even a sub-dominant mag- 
netic field is suscept ible to the heat-flux-driven buoyancy 
instab ility, or HBI (Quataert 2008 Parrish & Quataert 
20081. This instability, in the saturated phase, orients 
the magnetic field to be perpendicular to the tempera- 
ture gradient, a nd / C ond can be red uced to < 0.1 in the 
radial direction ( Parrish et al.|2009 ). On the other hand, 
the turbulent motion of gas can result i n more entan- 
gled mag netic field, restoring conduction ( Ruszkowski & 
Oh||2010|). For simplicity, we assume isotropic conduc- 



simulations with the suppression factor / con d = 0.1 and 
/cond = 0.3, which can be thought of as corresponding 
to the HBI dominant case and the turbulent case. To 
save computational time, these simulations are done us- 
ing 128 3 root cells (N root = 128). 

We find that thermal conduction does not significantly 
affect the overall evolution of the cluster, except that 
it slightly changes the temperature and density profiles 
when the cooling catastrophe occurs, and it can some- 
what change the timing of the cooling catastrophe itself 
by a few 10s of Myr. This result is consistent with what 
one would expect from Figure [6} the thermal conduction 
timescale i C ond is shorter than the cooling timescale at 
r > 100 kpc, where conduction can stabilize the local 
cooling instability, but at r < a few tens of kpc where 
the cooling catastrophe first starts to develop, i C ond is 
longer than i C ooi- 

To show the influence of thermal conduction near 



tion m our simulations, but we perform two additional 



r 100 kpc, we plot in Figure 11 the temperature pro- 
files from simulations with different suppression factors. 
When / C ond is larger, the thermal conduction is more ef- 
ficient, and thus the gas temperature is higher at r < 100 
kpc but slightly lower at r > 100 kpc because thermal 
energy is transported inwards from r > 100 kpc where 
the temperature peaks. Inside the core region, the tem- 
perature evolution looks very si milar with different f con d- 
Note that the outputs in Figure [TT] are shown 2 Myr af- 
ter the cooling catastrophe occurs in order to make a fair 
comparison in the central regions. While the / CO nd = 
and /cond = 0.1 runs collapse at almost exactly the same 
time (t — 295 Myr), the high conduction run actually col- 
lapses slightly earlier (about 20 Myr) because it damps 
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the initial temperature perturbations caused by the re- 
laxation at the beginning of the simulation which, in 
turn, are due to the slight deviation from hydrostatic 
equilibrium mentioned earlier. 

4.4. Supernovae la Heating 

The cooling catastrophe in these simulations takes 
place very close to the center of the BCG, where heating 
from Type la SN produced by stars in the BC G may have 
an effect on the formation of cooling flows (e.g. Domainko 
et al.|[2004| and can be important especial ly in the cen- 
tral regions (e.g. Conroy & Ostriker 2008). We do not 
include supernovae heating in our simulations but here 
we try to analytically estimate their contribution. Type 
II supernovae can be ignored in the BCG due to its rela- 
tively old stellar population. The SNIa ra te is ~ 0.1 per 



10 10 Mp) per century in galaxy clusters (jSharon et al. 



20071. Assuming each SNIa releases 10 ergs thermal 



energy into heating up the ICM, the heating rate then is 
expressed as M st eiiar (in M ) xlO 38 ergs yr _1 Mg 1 and 
the SNIa heating timescale can be estimated as 



iSNIa = 



P, 



gas 



y. 3 fcfoT 

2 ftrriH 



Psteiiar x 10 38 ergs yr 



-iM- 1 



(14) 



where p gas and p s teiiar are the local gas and the stellar 
densities. 

This timescale is plotted in Figure [6j which shows that 
although SNIa can be an important source of heating at 
r < 1 kpc and might slightly delay the formation of the 
cooling flow, it would not prevent the cooling catastrophe 
because t coo i is always shorter than isNia at all radii. 

4.5. Resolution and Numerical Method 

To test the influence of resolution on our results, we 
perform simulations with N mot = 64 and 128 (keeping 
a = —1.2) to compare with our standard -/V r0 ot = 256 
run. Figure 12 shows the density and temperature pro- 
files of the cooling gas at a similar stage of the evolu- 
tion (~ 2 Myr after the cooling catastrophe starts) . The 
overall behavior of the solution is seen to be relatively 
independent of resolution; however, there is a systematic 
trend for the temperature plateau in the intermediate 
region (between 100 pc and 20 kpc) to be even flatter 
as the resolution increases, while the transition radius 
shrinks. We do keep the maximum refinement level con- 
stant (/ max = 15), so the highest resolution achieved in- 
creases by a factor of two in each case; however, even 
for the lowest resolution run, this corresponds to 7 pc, 
considerably smaller than the ~ 100 pc transition radius. 
Therefore, we argue that it is not the maximum resolu- 
tion that is the key, but the resolution achieved in the 
early stages of the collapse. 

Since the point where the cooling catastrophe happens 
(the transition radius) is sensitive to the temperature 
and density profiles at earlier stages, we want to well re- 
solve the early evolution of the cooling gas in the center 
of the cluster. Changing a in the baryon mass refine- 
ment criterion (see Section 2.2) can affect the resolution, 



especially in the early stage inside the cluster core, and 
has a significant impact on the final results. The bottom 
two panels of Figure 12 show a comparison of runs with 
the same number of root cells iV root = 128, but different 



values for a. Again, the results are shown at a point 
roughly 2 Myr after the cooling catastrophe starts. At 
fixed radius, the gas density is lower and the tempera- 
ture is higher with more negative a, and the transition 
radius is smaller. Noticeably, with a — —0.1, the initial 
cooling region has a size of ~ 1 kpc, more than an order 
of magnitude larger than that in the run with a = —1.2. 
This is despite the fact that each calculation has the same 
maximum resolution. 

We also notice that there is a degeneracy between a 
and N mot : a higher N Ioot has a similar effect as a more 
negative a. This is because they both result in bet- 
ter mass resolution (see Equation [5]) in the center of 
the cluster as the cooling catastrophe develops. If we 
only increase the maximum refinement level l m ax with- 
out changing N root or a, the result does not change de- 
spite a smaller cell size at the highest refinement level. 
Thus we argue that it is crucial to have sufficiently high 
resolution at the early stage of the gas evolution. 

Although we have not been able to achieve complete 
convergence in these calculations, it is clear that higher 
resolution tends to produce flatter temperature profiles 
and smaller transition radii. 

Another difference between simulations with differ- 
ent resolution is that in low resolution runs (e.g. with 
N r oot = 64 and a — —0.2), the pressure drops dramati- 
cally inside r < 1 kpc when the central gas density and 
temperature show a sudden change, forming a pressure 
hole which then grows deeper and larger with time. The 
pressure hole is deeper when the resolution is lower. The 
gas inflow velocity is larger than that in the high resolu- 
tion runs and becomes supersonic at r ~ 1 kpc where the 
pressure gradient is the steepest, forming a sonic point 
and leading to a cooling catastrophe. Cold gas does not 
become rotationally supported as in the high resolutions 
runs but fragments inside the pressure hole via local cool- 
ing instability. 

Finally, we examine the impact of different methods 
for solving the hydrodynamics equations. As noted ear- 
lier, we use the Zeus method for our simulations because 
of its fast performance, and robust treatment of cold re- 
gions. To test the accuracy of its results, we also car- 
ried out test runs using the piecewise parabolic method 
(PPM) , which is third-order accurate with its high-order 
spacial interpolation. We used N roo t — 64, a = —1.2 
and found results which were in good agreement with 
the Zeus solver. 

4.6. Impact of the Choice of Initial Conditions 

In this section, we briefly examine how our results 
change as we vary details in the initial conditions. We 
start our simulations assuming the gas is in hydrostatic 
equilibrium. However, the NFW dark matter potential 
we use is fitted to the observed gas properties in the in- 
ner few hundred kpc region excluding the central ~ 10 
kpc where the BCG starts to dominate the potential, 
and thus the resulting NFW parameters can vary slightly 
depending on the boundary radius we choose when fit- 
ting the parameters. Therefore the gas is not in perfect 
equilibrium when the simulation starts, especially in the 
outskirts and inside < lOkpc. Examination of the early 
evolution of the cluster shows that the temperature and 
density of the gas in the central kpc increases by about 
30% and 50% respectively after about 20 Myr (~ t cross 
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Fig. 12. — The top two panels show the gas density (left) and temperature (right) approximately 2 Myr after the cooling catastrophe for 
three runs with varying root grid sizes, ranging from N roo t = 64 t o N roo t = 256 (but each run with a = —1.2). The bottom two panels 
show runs with N roo t = 128 but different a values (see Section|4.5|l, again for the same amount of time after the cooling catastrophe. 



at a few kpc) . They reach these values quickly, and then 
show little evolution until cooling becomes important on 
a few hundred Myr timescale. In order to check if the 
initial transients are important for the cooling, we also 
carried out a simulation where we let the system settle 
down for ~ 600 Myr (~ idyn at r ~ 100 kpc) without 
turning on cooling, so that the gas in the region of in- 
terest is in hydrostatic equilibrium when cooling starts. 
The later evolution of the cluster is the same as the run 
that starts with the cooling on. This indicates the tran- 
sients are not playing a role in cooling; however, it does 
mean that the true "initial state" of the core region (i.e. 
the density and temperature profiles after the transients 
die out) is slightly different from our given initial profiles. 
This is inevitable in the sense that the initial (observed) 



profiles are not in hydrostatic equilibrium. To test the 
general robustness of our results, we have also experi- 
mented with slightly different sets of NFW parameters 
for the dark matter and find that the gas properties at 
the outskirts are slightly different with different NFW 
parameters, but the cooling flow evolution in the core 
region is not affected. 

In many galaxy clusters, there is an offset between the 
X-ray emission center and the BCG, al though the offset 
tends to be smaller in CC clusters (e.g. | Sanderson et al. 
2009 1 . To see if the offset would significantly change the 



results, we have performed one simulation with an ini- 
tial offset of 20 kpc between the center of the gas and 
the gravitational potential. We find that the cluster gas 
settles down and re-centers on the BCG before the cool- 
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ing catastrophe happens. This is consistent with the fact 
that i CO oi > £dyn initially, and the cluster relaxes before 
cooling starts, and therefore the results do not differ from 
the simulations without the offset. 

To test the effect of other changes in our initial con- 
ditions, we performed one simulation without the initial 
random velocities. We find that changing the initial ran- 
dom velocity does not have a significant impact on the 
evolution of the cool core because the initial random ve- 
locity is damped before the cooling catastrophe happens. 

Since velocity perturbations do not directly perturb 
gas entropy, to further confirm that small scale pertur- 
bations do not grow outside the transition radius in our 
simulation, we performed a run with initial density per- 
turbations instead of velocity perturbations. To do this, 
we multiplied the density in each cell in the initial con- 
ditions by a Gaussian factor with a mean of unity and 
a standard deviation of 10%. Again, we do not see the 
grow th of any local instabilitie s . Th is is in agreement 
with Joung, Bryan & Putman (2011), who found that 
perturbations in a hydrostatic atmosphere did not cool 
unless the perturbation was sufficiently non-linear that 
the cooling time in the perturbation dropped below the 
time for the clump to accelerate to the local sound speed 
(roughly the dynamical time). 

We also performed a simulation without initial rota- 
tion. The gas in the very center in this run still eventu- 
ally becomes rotationally supported because the random 
initial velocities eventually are amplified due to the con- 
servation of angular momentum and a gas disk forms. In 
fact, even in our standard run, a smaller disk along the 
x-axis forms inside the major disk along z-axis, which can 
be seen in Figure [9j The size of the disk in the run with- 
out initial rotation is smaller at early times, and when 
the cooling catastrophe first occurs in that run, the gas 
in the very center has yet to become rotationally sup- 
ported; however, the required inflow velocity to balance 
cooling has already exceeded the sound speed, and so the 
flow passes through a sonic point (in the run with initial 
rotation, rotational support occurs before a sonic point 
develops). 

As noted earlier, the gravitational potential does play 
an important role in the cooling catastrophe and runs 
without a BCG produced significantly different results 
(see section 3.2 for more details). 

Finally, we carry out one simulation with a non CC 
configuration, where we use the same NFW dark matter 
profile for Perseus but set the initial temperature to be 
isothermal and compute the initial gas density assuming 
hydrostatic equilibrium. The initial t coo \ in the center is 
about 2 Gyr. Our simulation shows that after about 2 
Gyr, the cooling starts to run away and a cooling flow 
develops in a way which is quite similar to what we see 
in the simulations with initial CC configurations. The 
temperature plateau is slightly different, which we think 
has to do with the difference in the initial gas to dark 
matter ratio. We will return to this point in a later paper. 

4.7. Comparison with Previous Work 

In this section, we try to place these results in context, 
first making a link to steady-state cooling flow solutions, 
and then comparing to other (primarily simulation) work 
which looked specifically at only the developing cooling 
flow, and did not include feedback. 



The classic cooling flow model (Fabian 1994) predicts 
a "cooling flow" of 100s M /yr for rich clusters, assum- 
ing that in a steady state, without other heating sources, 
the gas flows inwards at a constant rate to replace the 
central gas that has cooled down and formed stars. The 
mass drop-out occurs over the central cooling-flow re- 
gion, and was originally assumed to cool and condense 
out in small clumps via a local cooling instability. This 
picture has been known to be in disagreement with ob- 
servations which usually indicate a star formation rate 
at least an order of magnitude lowe r than predicted b 
the steady state cooling flow model (Tamura et al.||2001 
O'Dea et al.]|2008l |Rafferty||McNamara fc Nulsen 2008| ) 
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Our simulations show that a steady state is not reached 
before the AGN feedback is potentially strong enough to 
balance cooling, and therefore we argue that this solu- 
tion is not relevant. However, it is interesting to see if 
we recover the steady-state result if we run the simula- 
tion for long enough. In Figure [13] we plot the gas inflow 
for a simulation with iV roo t — 64 that runs much further 
than our major runs. We find that after less than a hun- 
dred Myr, without any heating mechanism, the system 
approaches a steady state with roughly constant mass 
flow of M ~ 300 M Q /yr in the cluster core (r < 100 
kpc), consistent with the classic cooling flow prediction. 

Most previous simulation work on cool core clusters has 
focused on the heating process, especially AGN feedback, 
but they do usually include a pure co oling flow simulation 
in which heating is turned off (e.g. Croton et al. 2006). 
Our results are consistent with the se results inasmuch 
as there i s overlap. For example, Brighcnti & Math- 
( 2006 ) examined two dimensional models which also 
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found that a cooling-only model results in a relatively flat 
temperature profile (falling only by a factor of 2-3 over a 
range of 100 in radius). However, these simulations did 
not have the resolution (~ 1 kpc) to follow the cooling 
catastrophe in detail, and instead used a parameterized 
mass drop-out term in the mass-conservation equation. 
Similar results were found for a one -dimensional cooling 
flow in |Mathews fc Brighenti| ( |2003[ ). 

Our results are als o consistent w i th pre vious theoreti- 
cal work, for example, Bertschinger ( 1989 ) presented one- 
dimensional steady accretion models normalized to self- 
similar cooling waves that demonstrated slowly declining 
- or even increasing - temperature profiles. They even 
found a sonic radius (similar in radius to our transition 
radius) which occurred very close to the SMBH. 

Simulations which do not include a specific drop-out 
term and did not have high-resolution resulted in a much 
larger cooling region (a larger transition radius in our 
terminology) and a larger cooling flow immediately after 
t he cooling catastrophe t han found in our simulations. 



Kaiser & Binney ( 2003 ) provides a semi-analytic model 
of' the time evolution of the cooling flow assuming hydro- 
static equilibrium at every step of the evolution. Their 
density and temperature evolution agrees well with what 
we find for r > 1 kpc, when the inflow is subsonic and 
the gas is roughly in hydrostatic equilibrium. However, 
a semi-analytic model cannot describe what happens in- 
si de r < 100 pc wh ere cooling runs away. 



Guo et al. (2008) analyses the local and global insta- 
bilities in CC' clusters and shows that AGN feedback can 
suppress global radial instabilities. An interesting set of 
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Fig. 13. — The gas inflow of a simulation that goes ~ 100 Myr further than our standard run, reaching a steady state consistent with the 
classic cooling flow picture. The left panel shows the gas inflow profile at each output time and the right panel shows the time evolution of 
the inflow at a radius of 3.42 kpc. The radius is c hosen to be large enough to avoid the fluctuation in the very center, while still providing 
a good estimate of the mass flow (see Section |4.7[l . 



recen t simulations (McCourt et al. 2011 Sharma et al. 
20111 have addressed the roie of local cooling instabili 



ties m simulations of a cool core, suppressing the cooling 
catastrophe by instituting a global heating term. They 
showed that local cooling instabilities can form as long 
as the timescale for a cooling instability is less than 10 
times the dynamical time. That result does not conflict 
with our finding in this paper that local cooling instabil- 
ities do not occur before the global instability (or catas- 
trophe) because of the global heating term instituted in 
those paper. 

Finally, we note that we do not include feedback and 



model AGN energy injection 



ian et al. | |2010p , althoug 



e.g., 


Omma et al.| 


2004 


Dubois et al. 2010| Fab- 



work. Also, we have focused on the evolution of cool core 
clusters rather than the formation of cool core clusters, 
and so it remains possible that isotropic thermal conduc- 
tion can play a role in th e earlier evolution of a non-cool 
core clusters ( Voit 2011). 



5. CONCLUSION 

We performed high resolution simulations of an ideal- 
ized, cool-core X-ray cluster loosely based on Perseus, in 
order to better understand how gas cools and accretes on 
to the SMBH. The use of an AMR code was essential in 
order to resolve the important scales, ranging from the 
Mpc virial radius of the cluster, down to the pc scale, 
which is the size of the cold gas clumps. We do not in- 
clude feedback from the AGN and so the current work is 
not a complete description of cool-core cluster evolution; 
however, it does shed light on where cold gas condenses 
out of the flow, where heating needs to occur, and what 
observational signatures we would predict for such an 
object. 

We simulate an idealized cluster, but include small- 



scale velocity perturbations, as well as some net rota- 
tion, as implied by cosmological simulations. We include 
radiative cooling, but no heating term. The simulation 
resolves regions of high density and short cooling time. 
We present our primary results below. 

• We find that as the cluster core cools, the temper- 
ature profile remains remarkably fiat, from a few 
hundred pc out to at least 20 kpc, with a change 
of only a factor of ~ 2 in temperature over this 
range. We argue that this occurs because the gas 
flow adjusts such that adiabatic compression bal- 
ances cooling at each radius over this range. The 
gas that does cool does so via a global cooling catas- 
trophe that occurs first at a 'transition' radius of 
about 10 pc from the SMBH in our best resolved 
simulation. This occurs about 300 Myr after the 
start of the simulation, which is about the cooling 
time of the gas in the initial configuration. 

• We stop the simulation about 20 Myr after the 
cooling catastrophe when the accretion rate onto 
the black hole (which we measure at a radius of 
400 pc) becomes sufficiently large that a reasonable 
feedback efficiency could balance radiative cooling 
in the flow. At this point, when we examine the 
amount of cold gas, there is a distinct lack of gas 
below a few keV (because of the flat temperature 
profile mentioned above) , and is in agreement with 
observational constraints on Perseus. 

• Although the simulation is intrinsically three- 
dimensional, the solution outside of the transition 
radius is well-described by spherically symmetric 
flow. In particular, we find that local thermal 
instabilities do not grow in the cool-core cluster 
within 300 Myr (i.e. before the global cooling 
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catastrophe) , at least outside of the transition ra- 
dius around ~ 100 pc. This is despite the fact that 
we seed the initial flow with substantial perturba- 
tions. 

• Isotropic thermal heat conduction does not signifi- 
cantly affect this result. It does not lead to a large 
delay, nor does it significantly affect the density and 
temperature profile in the center. Heating from 
Type la SN is similarly unimportant. 

• The final result is not very sensitive to the gas ini- 
tial conditions, but is sensitive to the presence of 
the BCG. Without a BCG, the gas temperature 
drops much more rapidly with radius, in conflict 
with observations. 

• We find that some aspects of the solution are quite 
sensitive to the resolution of the simulation. In 
particular, we find that very high mass resolution 
is required in the gas near the cooling flow, and 
that lower resolution runs resulted in both a steeper 
slope in the temperature profile in the intermediate 
range from the transition radius to 20 kpc, as well 
a larger initial transition radius. 

These results are intriguing for a number of reasons. 
First, they imply that the real "cooling-flow" problem 
arises very close to the center of the cluster (i.e. within 
a few hundred pc), in the sense that it is only in this 
region that gas cools to temperatures below those ob- 
served. Therefore, the lack of cooler gas seen in X-ray 
spectroscopy may arise not because of the way that heat- 
ing operates, but because of the temperature plateau 
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caused (we argue) by the BCG potential well. Mass drop- 
out occurs only at the very center of the cluster. 

The fact that mass drop-out occurs first very close to 
the SMBH also provides a mechanism for thermal regu- 
lation to operate. One of the unanswered questions for 
AGN regulation is how the properties of the gas in the 
entire cooling region, which stretches out to 100 kpc, can 
control what happens in the center hundred pc. These 
simulations help to answer this by showing that this nat- 
urally occurs - the lack of local thermal instabilities in- 
sures that the flow is coherent, and, as we show, the 
result is cold gas condensing out only inside the tran- 
sition radius, where it can almost immediately form an 
accretion disk and feed the SMBH. 

Of course, this doesn't relieve us of the need for en- 
ergetic feedback, as we show that the mass inflow rate 
quickly grows and after a few hundred Myr, reaches hun- 
dreds of solar masses per year, in agreement with steady- 
state models, and in disagreement with observations. We 
argue that, before this occurs, the accretion disk that we 
see feeds gas to the SMBH, generating 'radio-mode' feed- 
back that heats up the gas. How exactly this feedback 
occurs is not clear and not addressed in this paper. We 
will examine this point in more detail in future work. 
Another aspect which is not explained by these results 
is the generation of Ha filaments seen in most cool core 
clusters. Here we simply note that they do not naturally 
occur at large radii via a local cooling instability before 
the flow becomes globally unstable. 
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